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LES, DNS and RANS for the Analysis of 
High-Speed Turbulent Reacting Flows 


P.J. Colucci, F.A. Jaberi and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York 
Buffalo, NY 14260-4400 


Abstract 

The purpose of this research is to continue our efforts in advancing the state of 
knowledge in large eddy simulation (LES), direct numerical simulation (DNS) and 
Reynolds averaged Navier Stokes (RANS) methods for the computational analysis of 
high-speed reacting turbulent flows. We have just finished the third year of Phase II of 
this research. This report provides a detailed description of our most recent findings 
in the research supported under this program. 


Technical Monitor: 


Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, Mail Stop 197, Tel: 
804-864-2298) is the Technical Monitor of this Grant. 
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1 Introduction 


We have just finished Year 3 of Phase II activities on this NASA LaRC sponsored project. 
The total time allotted for this phase is three years; this phase was followed at the conclusion 
of Phase I activities (also for three years). Thus, in total we have completed 6 years of LaRC 
supported research. A proposal for continuation of these efforts is pending at the NASA 
LaRC. 

Our work in this project can be divided into the following categories: (1) development 
of LES methodologies via PDF methods and numerical solution of the PDF via Monte 
Carlo schemes. (2) Development of algebraic turbulence closures for high speed reacting 
flows. (3) Implementation of the models developed in (2) in the computer code LAR.CK. 
(4) Applications of the LES-PDF methodology for flows of interest to NASA LaRC. (5) 
Utilization of DNS for model assessments. We do not have significantly new results since 
our last report pertaining to item (2) to present here. We have made significant progress in 
regard to tasks in item (3). However, since a thesis is under preparation to be submitted to 
NASA LaRC within the next month, this item is not discussed in the present report. Item 
(4) is relatively new and no results are presently available to present. The major thrust of 
this report is to provide a detailed descriptions of our achievements under items (1) and (5). 


2 The Filtered Density Function for LES of Chemi- 
cally Reacting Flows 

The prediction of turbulent flows has been the focus of extensive research due to its impor- 
tance in practical engineering applications. In past decades the analyst has had three main 
tools in the simulation of turbulent flows. The traditional approach in prediction of flows for 
engineering applications has been to consider the Reynolds averaged Navier-Stokes (RANS) 
equations. In RANS computation, the equations of motion for time averaged transport vari- 
ables are solved. As a consequence of the averaging procedure, models must be provided to 
account for the effects of the turbulence on the mean flow. Alternatively, direct numerical 
simulation (DNS), in which all scales of the flow are resolved, offers the advantage that the 
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turbulence is explicitly calculated rather than modeled. Due to the very fine grids required in 
such model free simulations, DNS is typically restricted to very simple flows making its appli- 
cation to engineering problems impractical. Large eddy simulation (LES) offers a compromise 
between DNS and RANS computation (Givi, 1989; Libby and Williams, 1994; Givi, 1994; 
Rogallo and Moin, 1984; Lumley, 1990; Galperin and Orszag, 1993; Voke and Collins, 1983; 
Reynolds, 1990; Oran and Boris, 1991). In LES, local volume averaging is utilized to ob- 
tain a smoothed solution which can be captured on reasonable grids. In contrast to RANS 
computation in which all time (and implicitly length) scales must be modeled, only the 
smallest scales of the turbulent flow require closure, while the larger turbulent flow structure 
is resolved. 

Over the past thirty years since the early work of Smagorinsky (1963) there has been rel- 
atively little effort, compared to that in RANS calculations, to make full use of LES for 
engineering applications. The most prominent model has been the Smagorinsky eddy vis- 
cosity based closure which relates the unknown subgrid scale (SGS) Reynolds stresses to 
the local large scale rate of flow strain (Lilly, 1965; Lilly, 1967). This viscosity is aimed to 
provide the role of mimicking the dissipative behavior of the unresolved small scales. The 
extensions to “dynamic” models (Germano, 1992; Germano et al , 1991; Moin et a/., 1991; 
Lilly, 1992) have shown some improvements. This is particularly the case in transitional 
flow simulations where the dynamic evolutions of the empirical model “constant” result in 
(somewhat) better predictions of the large scale flow features. 

A survey of combustion literature reveals relatively little work in LES of chemically reacting 
turbulent flows (Givi, 1989; Pope, 1990). The primary stumbling block in this utilization 
is the lack of closures which can accurately account for the influence of unresolved subgrid 
scale (SGS) fluctuations (Madnia and Givi, 1993). It appears that Schumann (1989) was one 
of the first to conduct LES of a reacting flow. However, the assumption made in this work 
simply to neglect the contribution of the SGS scalar fluctuations to the filtered reaction rate 
is debatable. The importance of such fluctuations is well recognized in RANS computation 
of reacting flows in both combustion (Libby and Williams, 1980; Libby and Williams, 1994; 
Jones and Whitelaw, 1984; Jones, 1994) and chemical engineering (Brodkey, 1975; Toor, 
1975; Hill, 1976; Brodkey, 1981) problems. Therefore, it is natural to believe that these 
fluctuations will also have a significant influence in LES. 
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Modeling of scalar fluctuations in RANS has been the subject of intense investigations since 
the pioneering work of Toor (1962). The aim of statistical moment methods is to provide a 
closure for these correlations in terms of the mean flow variables. Eddy break up type closures 
(Spalding, 1977), which were popular in the 1970’s, have achieved only a limited degree of 
success. Because of the lack of models with universal applicability to accurately predict the 
scalar correlations in turbulent reactive flows, simulations involving turbulent combustion 
are often met with a degree of skepticism. Another approach which has proven particularly 
useful is based on the probability density function (PDF) or the joint PDF of scalar quantities 
(O’Brien, 1980; Pope, 1985; Dopazo, 1994). This approach offers the advantage that all the 
statistical information pertaining to the scalar field is embedded within the PDF. Because of 
this feature, PDF methods have been widely used in RANS for a variety of reacting systems 
(see Dopazo (1994) for a recent review). The systematic approach for determining the 
PDF is by means of solving the transport equation governing its evolution (Lundgren, 1967; 
Lundgren, 1972; O’Brien, 1980). In this equation the effects of chemical reaction appear 
in a closed form. However, modeling is needed to account for transport of the PDF in the 
composition domain of the random variables. In addition, there is an extra dimensionality 
associated with the composition domain which must be treated. An alternative approach 
is based on assumed PDF in which the PDF is parameterized a priori in terms of its lower 
(usually the first two) moments. Obviously, this method is ad hoc but it offers more flexibility 
than the first approach. Presently the use of assumed methods in RANS is justified in 
cases where there is strong evidence that the PDF adopts a particular distribution (Rhodes, 
1975; Jones and Priddin, 1978; Janicka and Kollmann, 1979; Lockwood and Moneib, 1980; 
Bockhorn, 1988; Priddin, 1991; Madnia et al., 1992; Frankel et ai, 1993b; Miller et at ., 1993). 

In an attempt to account for the effects of the subgrid scales on the filtered reaction rate, 
the eddy break up concept has been adopted for use in LES (Fureby and Lofstrom, 1994; 
Garrick, 1995). As in the case of RANS computation, such closures have achieved only a 
limited degree of success. Despite the demonstrated capabilities of PDF methods in RANS, 
their use in LES is limited (Givi, 1989; Madnia and Givi, 1993). The first application of 
PDF-LES is due to Madnia and Givi (Madnia and Givi, 1993) in which the Pearson family of 
PDF’s are used for modeling of the SGS reactant conversion rate in homogeneous flows under 
chemical equilibrium conditions. This procedure was also used later by Cook and Riley (Cook 
and Riley, 1994) for LES of a similar flow. The extension of assumed PDF model for LES 
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of nonequiUbrium reacting shear flows is reported by Frankel et al (Frankel et al, 1993a; 
Frankel, 1993). While the generated results are encouraging, they do reveal the need for 
more systematic schemes. Most of the drawbacks of these schemes can be overcome by 
considering the solution to the equation governing the evolution of the PDF. Because of the 
added dimensionality of the PDF transport equation, its solution by conventional numerical 
methods is possible in only the simplest of cases (Janicka et al, 1979). An analysis performed 
by Pope (1981) suggests that the solution of the joint velocity-scalar PDF by finite difference 
methods is impractical for more than three scalars. However, the Lagrangian “Monte Carlo” 
scheme (Pope, 1985) can be used for this purpose. While PDF transport methods have 
enjoyed significant development in the past decade, no attempt has ever been made in its 
utilization for LES. 

In this work, the Lagrangian PDF methodology is utilized as a means of providing the scalar 
SGS correlations for chemically reactive flows under nonequilibrium conditions. The La- 
grangian approach, which has proven extremely effective in Reynolds averaging procedures, 
offers an attractive means of predicting the evolution of the “filtered density function” (FDF), 
which is essentially the PDF of SGS variables. In this paper we demonstrate its enormous 
potential for LES of turbulent combustion. 


3 Mathematical Formulation 

3.1 Governing transport equations 

In the present treatment, incompressible flows undergoing isothermal reaction are considered. 
The restrictions regarding variable density due to exothermicity require attention, however 
they do not significantly alter the fundamental analysis and are readily removed. In the 
computational treatment of such flows, the primary independent transport variables are: 
the velocity vector tq (i = 1,2,3), the pressure p, and the species mass fractions <f> a (a = 
1,2, .. . , N s ). The conservation equations governing these variables are given by: 
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Continuity: 


( 1 ) 



Conservation of momentum: 

duj duiUj dp drij 

dt + dx{ dxj dx{ 


( 2 ) 


Conservation of chemical species: 

d<f> a dui<f> a 
dt dxi dxi 


( 3 ) 


The viscous stress tensor T tJ - and mass flux Jf of chemical species a are given by: 


Tij = P 




( 4 ) 


r d<t>* 

dx { ' 


( 5 ) 


where p is the molecular viscosity, T = p/Sc and Sc is the Schmidt number. The transport 
equations given above provide a complete description of the reactive system. However, the 
variations in length and time scales would require computational resolutions which are too 
prohibitive for even the fastest of today’s supercomputers(Zang et al., 1989). In RANS, solu- 
tions of the time averaged transport equations are attempted. In LES, instead of averaging 
over all time (and implicitly length) scales, the transport pertaining to the larger, energy 
containing eddies are considered. This allows resolution of the lower frequency turbulent 
structures. Operationally, this involves the use of a local “spatial filter” (Aldama, 1990): 


/ +oo 

f(x',t)G(x’ -x)dx' 

-CO 


( 6 ) 


where / represents the filtered value of the transport variable /; G(x) denotes the filter and 
f = f -J denotes the fluctuations from the filtered value. With the application of the 
filtering operation to the governing transport equations, we obtain: 


dip 

dxi 


= 0 


( 7 ) 
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( 8 ) 


duj duiUj 
dt dxi 


dp | dr l3 dr t] 

dx } dxi dxi 


d <k = dJ]_ m I z 

dt dxi dxi dxi 

where Tij = uju] - UiUj , M“ = - Ui<j> a and u> a are unclosed terms. 


( 9 ) 


3.2 Modeling of unresolved scales 

The quantities T \j and M“ introduced in the previous section represent the effects of the 
small scale turbulence on the larger scales. A model must be introduced to properly account 
for these effects. The use of gradient diffusion models, particularly the one introduced by 
Smagorinsky (1963), have enjoyed a reasonable amount of success. While more sophisticated 
two-equation models, Reynolds stress models and Algebraic Reynolds stress models have 
undergone significant development for RANS computation, the turbulence models developed 
for LES are considerably less complicated. This is at least partially due to the fact that in 
LES less of the turbulence is modeled (and more is resolved) than in RANS computation. 
The Smagorinsky model is a gradient diffusion model which takes on the form: 

Tn - (Si S /3)T kh = -2/ttSij (10) 


where Sij is the large scale strain rate tensor and 

fh = C.Aly/SijSij. ( 11 ) 

The dynamical interaction between resolved and SGS fields in isotropic turbulence was stud- 
ied in detail by Kerr et al. (1996). Their results at different Reynolds number show that 
there is a strong correlation between large scale energy and SGS dissipation. This correlation 
implies that SGS dissipation is most effective at those physical locations which are rich in 
the large-scale energy. They conclude that at least in isotropic turbulence large-scale energy 
is a good indicator of nonlinear SGS interactions. Based on these observations we propose 
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to model the eddy viscosity with the following (MKEV) model: 


(*t 


Ck&G\f\ 


UiUi — UiUi\ 


( 12 ) 


where the carat denotes the filter at the secondary level which has a characteristic size larger 
than that of grid level filter, and Aq is the filter width. We have found reasonable results for 
A g/Aq = 3 where Aq is the characteristic size of the secondary level filter. Note that the 
observations by Kerr et al were for homogeneous flows and the model has been presented 
so as to satisfy Galilean invariance for inhomogeneous flows. For homogeneous isotropic 
flows Hi — > 0 as A g oo. This model is essentially a modified version of that proposed by 
Bardina et al. (1983), which utilize equal sizes for the grid and secondary filters. 


A similar model may be used to close the subgrid mass fluxes: 


M" = -I\ 


dx{ 


where T* = p. t j Sc t , and Sc t is the turbulent Schmidt number. 


(13) 


Thus far we have not yet addressed the issue of how to deal with scalar correlations in 
the filtered chemical source terms. For second order isothermal reactions, the correlation 
4> a <t>p must be properly modeled in order to correctly account for its effects on the larger 
scales. The SGS unmixedness (r a p = <j>ot<t>p — <f> Q <f>p) represents the effect of SGS scalar-scalar 
corelation on the filtered chemical source terms. In RANS calculations the importance of 
the unmixedness term has long been recognized. In LES, the importance of this term has 
not been well appraised. While the SGS stress and flux terms discussed in this section are 
of a convective nature and can be reasonably well modeled by a diffusive process, the same 
cannot be said for the unclosed terms in the species production rates. Because the physical 
mechanism of the SGS stresses and fluxes is inherently different from the scalar correlations 
in the chemical source terms, it is expected that the models will differ. In fact, when 
eddy viscosity concepts are extended to treat chemical source terms, the resulting models 
(“eddy break up models”) perform mediocre at best. The focus of the following sections 
is to discuss how the LES-PDF methodology is used to overcome the closure problem of 
the chemical source terms and to develop robust numerical methods for the simulation of 
turbulent reactive flows. 
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3.3 The filtered density function (FDF) 


The most effective means of dealing with the statistical closure problem is to consider the 
joint PDF of the transport variables. For the purpose of simplicity, here we first consider the 
scalar array <p(x,t) only; treatment of joint velocity-scalars follow a similar procedure. For 
this array, the “filtered density function” (FDF), denoted by Pl, is defined as (Pope, 1990): 

/ +oo 

Q [V’) < / > (x / , t)j G(x' — x)dx' (14) 

* oo 

N 3 

g[lp,(p(x,t)} = Sty - 0(x,*)] = H 5[lp a - 4>a{x,t)} (15) 

a=l 

where 8 denotes the delta function and ip denotes the “composition domain” of the scalar 
array. The term g[<p — ip(x,t)} is the “fine-grained” density (O’Brien, 1980; Pope, 1985), and 
Eq. (15) implies that the FDF is the spatially filtered value of the fine-grained density. The 
FDF in this way, was first defined formally by Pope (Pope, 1990) and it has been discussed 
by others (Madnia and Givi, 1993; Gao and O’Brien, 1993). However, it has yet to be 
systematically used in actual LES. 

Evaluations of spatial filtered values of the transport variable are achieved by integrating 
the FDF, exactly similar to that in PDF. For example, for the function A(<p), the filtered 
(desired LES) value is: 

/ +oo f+oo 

A{x',t)G(x' -x)dx' = / A{ip)P L (ip]x,t)dip. (16) 

-oo — oo 

More generally, the conditionally filtered value of any function A(x, t) may be evaluated by 
integration over the whole of the compositional domain: 

/ +oo 

-4(x', 0 - - x)dx'. (17) 

-OO 

The unconditional average may be obtained by: 

/ +oo 

(A(x,t)\ip) P L (ip-x,t)dip. (18) 

-OO 
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If the variable A is completely determined by the compositional variable ip, then this last 
expression reduces to Eq. (16). 


3.4 FDF transport equation for reactive scalars 


The equation governing the transport of the scalar PDF may be derived by averaging the 
evolution equation of the fine-grained PDF (Lundgren, 1967; Lundgren, 1972; Stratonovich, 
1963; Pope, 1976). In an analogous manner, the FDF transport equation may be arrived 
at by spatially filtering the same fine-grained transport equation. Alternately, Pope (1985) 
suggested an approach in which two independent expressions for ( DQ/Dt ) are equated. The 
same methodology may be applied to derive the evolution equation of the FDF. Another 
procedure, based upon Eqs. (14) and (17), is presented here. 


Differentiation of Eq. (14) with respect to time yields: 

dP L {ip-,x,t) [<*> d<p a d8((p(x' ,t) - ip) 


dt !-oo dt 

d f°° d<j>, 


dfio 


G(x' — x)dx' 


/ -~6{(p{x! ,t) - ip)G(x' -x)dx'. 


dlpa J-oo dt 

This expression can be combined with Eq. (17) to obtain: 


dP L (ip\x,t) d 


dt 


dlpa 


' d(p a 
, dt 


\xp) P L (ip-,x,t) 


Utilization of Eq. (3) into this last expression yields: 


dP L {ip;x,t) d 




dt dlpa 

The convective term may be decomposed into large and small scale contributions: 


d 


dUi<p 0 

dip a [\ dxi 


\ip) P L {ip-,x,t) 


d 


d*Pa 




Ui ^SL\A -u,( d -^-\d 


dxi 


dxi 


P L (tp;x,t) 


(19) 


( 20 ) 


P L (ip-,x,t )\ . (21) 


(22) 
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A similar manipulation to that used to obtain Eq. (19) can be used to recast the large scale 
convective term in a more familiar form: 


3 


Ui 



P L (^;x,0 


duiP L (V>;x,t) 

dx{ 


(23) 


The FDF transport equation then takes on the form: 


DPl duiPi 
dt dx{ 



d[o; a (V>)P L ] 

dtp a 

(24) 


The first two terms on the left hand side of the equation represent the large scale convection 
of the FDF in physical space and appear in closed form. The first term on the right hand 
side denotes the effects of unresolved small scale convection of the FDF in physical space 
and requires closure. The second term on the right hand side represents molecular mixing 
and diffusion and requires modeling. In general, the mixing term tends to homogenize the 
fluid and hence lowers the scalar variance. The last term on the right hand side is due to 
chemical reaction which is in a closed form. The major advantage of PDF methods is that 
the scalar correlations in the filtered reaction rate do not need to be modeled. Instead, 
the effects of reaction are exactly accounted for without approximation. Note that the last 
two terms contain derivatives in compositional space rather than physical space. This is 
reflected by the fact that the processes of mixing and chemical reaction serve to change the 
compositional makeup of the mixture rather than to provide a mechanism for transport in 
physical space. Effectively, the last two bracketed terms on the right hand side are fluxes in 
the compositional domain. 


3.5 Modeled FDF transport equation 

While no modeling is required to account for the reaction source term in Eq. (24), the small 
scale convective term must be modeled to account for the effects of the turbulent transport 
of the FDF in the physical domain, and a model must be presented for the conditional 
expected diffusion in the molecular mixing term to account for transport in the domain of 
the compositional variables. 
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A gradient diffusion model of the form 


d 

di>a 


d<j> ( 

u,-x~ 



- Ui 


d_K 

dxi 




dxi 1 dxi 


(25) 


is adopted to account for the small scale convective effects. This equation may be integrated 
to generate moment expression for the subgrid moments. For example, the first subgrid 
moment of Eq. (25) is given by: 


die _ d<f>, 

Ui~ Ui 


dxi 


dxi 


d r 

dxi t dxi 


(26) 


Similar expressions may be derived for the higher order subgrid moments. One must be 
careful in the integration procedure to obtain Eq. (26) in the case where the spatial filter 
varies with respect to the spatial coordinate. This is due to the fact that for nonuniform 
filters, the filter and derivative operators do not commute. 


Modeling of the molecular mixing term has been the focus of intense investigation (Pope, 
1985). Most of the models currently in use are based on the family of Coalescence/Dispersion 
(C/D) closures (Curl, 1963; Dopazo and O’Brien, 1976; Pope, 1976; Janicka et a/., 1979; 
Pope, 1979; Pope, 1982; Kosaly and Givi, 1987; McMurtry and Givi, 1989). Our previous 
investigations (Kosaly and Givi, 1987; McMurtry and Givi, 1989) indicate that in homo- 
geneous turbulence the predicted results are sensitive to the choice of the mixing model. 
However, it is anticipated that in the context of inhomogeneous LES, different members of 
C/D models behave very similarly, particularly for the lower order moments. Amongst these 
members the one which is the most convenient to use is the linear mean square estimation 
(LMSE) model (O’Brien, 1980; Dopazo and O’Brien, 1976). In the adaptation of this model 
for LES, the value of the transport quantity relaxes to the “spatially filtered” value at the 
point of interest. For a scalar governed by Fickian diffusion, the molecular mixing within 
the subgrid is modeled as: 


d 




0 ( r W\ + a 


ra.fi 


dxi \ dxi I dip a [ 2 




(V>o 


4>*)Pl 


(27) 


where Cl m is the frequency of mixing within the subgrid and can be related to the subgrid 
diffusion coefficient and the filter length: fl m — Co(r + r t )/A^. 
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With the closures given by Eqs. (25) and (27), the final modeled FDF transport equation is 
given by: 


dP L d[ujP L ] 
dt + dx t 


±< r + T ^ + — 
dxi * dxi di/ic 


CM 


05 ‘m 


(0. - MPl] ~ 


d\uJ^p)P^ 

dV’a 


(28) 


3.6 Consistency between FDF and moment closure approaches 

The modeled evolution equation for the FDF, Eq. (28), may be integrated to obtain transport 
equations for the moments. While the derivations to this point have been completely general 
with regard to variable spatial filters, the restriction of uniform filters is invoked for this 
section. With a few notable exceptions, nearly all SGS modeling has been developed for 
uniform filters. 

In the case that the spatial filter is not a function of space the convective effects may be 
taken into account via the decomposition: 

(Ui IVO Pl = UiP L + [(Ui\ij>) - Ui]P L = u{Q + [uig — Ui~g] . (29) 

This result is easily derived by imposing the restriction of constant filter size on the decom- 
position given by Eq. (22). The attractiveness of the decomposition given by Eq. (29) is 
that upon integration, it yields results identical to those in conventional LES. For example, 
the first moment of Eq. (29) is: 


Ui<f> Q = Ui<t> a + [Ui<j) a - Ui<j> a \. 


(30) 


The term in brackets in Eq. (30) is the generalized scalar flux in the format defined by 
Germano (1992). Thus, if a gradient diffusion model of the form 


l(w.lV’) 


U; 


i]P L = -r, 


OPl 

dxi 


(31) 


is adopted, the first moment of the FDF satisfies: 


\u. 


<f> a - Ui(f> a ] = -Tt 


dxi 


(32) 
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This is similar to the Smagorinsky type closure (Lilly, 1965; Lilly, 1967) as used in conven- 
tional LES (with or without the dynamic (Germano, 1992; Germano et a/., 1991) model). 
Of course, this similarity does not imply that the FDF is “equivalent” to the Smagorinsky 
closure. Although the two closures yield similar results for the “first subgrid” moment of 
a “non-reacting” scalar, the FDF includes all the higher moments of all the reacting scalar 
variables. 


For uniform filters, the conditional expected diffusion can be decomposed into diffusive and 
dissipative parts: 


d 
dip o 




s_ ( dPL\ _ y 

dxi l dxi I dtp a dlp 0 


'^d<p a d(p Q[ \ ' 


(33) 


This isolates the effect of scalar dissipation for which LMSE provides a closure: 


d 2 


dfpadtpo, 


\,dp a d(f>c t \ 


d 

dip a 


cm 




(*> - K) ft] ■ (34) 


This again, yields results consistent with conventional LES for the first two subgrid moments. 
While the first moment of Eq. (34) yields a trivial result (molecular mixing has no effect on 
the mean), its second moment generates an expression for the dissipation: 


dxi dx t 


— — 2 . 

o Kvi - <Pa )• 


(35) 


The modeled FDF transport equation, Eq. (28), may be integrated to obtain transport 
equations for all the desired moments. For instance, with the assumption of a uniform filter 
width, the equation for the first subgrid moment is given by: 


du -jfa 
dt dx t 


^( r+r *> 


dxi 


+ 


Similarly, the transport equations for higher order moments may be derived: 


(36) 


d°l dujal 
dt dxi 


d_ 

dx t 


< r+r '>S 


C^mCr 2 T 2(r + 


U) 


dpa d(pot 
dx, dxi 


d - 2(^> a cj Q 


<t> a u*) (37) 


where cr 2 = (p^ — <p Q is the generalized subgrid variance of chemical species <j> a . It should 
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be noted that the quantity a \ is more useful in LES than the quantity <t>' a <)>' a as the former 
term represents the entire deviation between the subgrid and resolved scales while the latter 
term does not account for cross and Leonard stresses. Under the restrictions of no chemical 
reaction, it is clear that direct solution of the moment equations (36) - (37) yields the 
exact same information for the first two moments as those generated by the solution of the 
modeled FDF transport equation. This consistency is exploited by comparing the solution to 
the moment equations and the FDF under non-reacting conditions. This allows verification 
that the numerical solution to the FDF is accurately capturing the effects of convection, 
mixing and diffusion. 

While it has been shown that consistent models for turbulent transport and molecular mixing 
(i.e. dissipation) may be derived, no such consistency condition exists for the reaction source 
term. This is a consequence of the simple fact that no model for reaction is required in the 
FDF treatment, and hence the analogous moment closure does not exist. Under chemically 
reactive conditions, the advantage of the FDF approach is obviated; no model is required to 
treat the nonlinear (and arbitrarily complicated) chemical source terms, while the moment 
equations require such a model. This is particularly important as at the present time no 
generally applicable reaction model has ever been proposed for use in LES. 

3.7 The Lagrangian approach, Langevin equation and equivalent 
systems 

The solution of Eq. (24) provides all the statistical information pertaining to the scalar 
array </>(x,f). This equation can be solved most effectively via the Monte Carlo scheme 
(Pope, 1981). Currently, two classes of Monte Carlo schemes are available: Eulerian and 
Lagrangian. In the Eulerian scheme, the PDF within the subgrid is represented by an 
ensemble of computational elements (or particles) at “fixed” grid points. These elements are 
transported in the “physical space” by the combined actions of large scale convection and 
diffusion (molecular and subgrid). In addition, transport in the “composition space” occurs 
due to chemical reaction and subgrid molecular mixing. In previous unpublished work, the 
Eulerian Monte Carlo method was tested for LES of a non-reacting shear flow. Expectedly, 
the results were not encouraging. The major difficulty with the Eulerian formulation lies in 
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the numerical implementation of the large scale convection. Due to “grid-based” nature of 
the scheme, excessive artificial diffusion is created which degrades the solution very notably. 
It is crucial to realize that the errors induced by this scheme are not due to the FDF 
formulation, but rather to the numerical implementation of the large scale convection. 

A remedy for the problem noted above is to divorce from the Eulerian discretization and to 
invoke the Monte Carlo solver in a “grid free” Lagrangian manner. The Lagrangian procedure 
involves transport of the elements within the “whole” computational domain of interest. 
The advantages of this procedure in reducing the amount of numerical diffusion in DNS 
are well- recognized (Chorin, 1968; Chorin and Marsden, 1979; Leonard, 1980; Majda, 1988; 
Sarpkaya, 1989; Ghoniem, 1991; Gustafson and Sethian, 1991). The basis of the Lagrangian 
solution of the FDF transport equation relies upon the principle of equivalent systems (Pope, 
1985; Pope, 1994). Two systems with different instantaneous behaviors may have identical 
statistics and satisfy the same FDF transport equation. In the Lagrangian solution procedure 
each of the particles obeys certain equations which govern its transport. It is important to 
recognize that these particles are not fluid elements. In fact, while fluid parcels follow smooth 
trajectories, it will be shown that the Monte Carlo particles follow trajectories which are 
continuous but not differentiable. The significance of these notional particles is that they 
are developed in such a way that they evolve with the same collective statistics as genuine 
fluid particles. 

The Monte Carlo particles undergo motion in the physical space by convection due to the 
filtered mean flow velocity, and diffusion due to molecular and subgrid diffusivities. The 
general diffusion process may be represented in a stochastic manner by the Langevin equation 
(Pope, 1985; Risken, 1989; Gardiner, 1990; Gillespie, 1992): 

dXi(t) = A(X(f), t)dt + t)dWi{t) (38) 

where Xi is the Lagrangian position of a stochastic particle. The entities D, and B are known 
as the “drift” and “diffusion” coefficients, respectively, and W{ denotes the stochastic Wiener- 
Levy (Karlin and Taylor, 1981) process. This process is best understood by considering the 
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function Wd{t n ) which changes at discrete time intervals: 

Wrf(<n) = (Ai) 1/2 i;6, (39) 

t=i 

where £ n (n = 1,2, ...,N) are N independent normalized Gaussian random variables and 
the time interval [0, T ] is divided into N equal subintervals of duration At = T/N. Consider 
the increment 

AW d (t n _ 1 ) = U&t) l/2 - (40) 

The Weiner process can be defined as Eqs. (39) - (40) in the limit At — * 0. Note that 
although the process is continuous, it is not differentiable since AWj,/ At is undefined as At 
vanishes. For the system considered here, the stochastic particle locations are governed by 
the solution to the following stochastic differential equation: 


dXi(t) 


Ui + 


3(r + r,) 


dxi 


dt + wr + Tttf^dWi. 


(41) 


It must be emphasized that although the Langevin equation given by Eq. (41) is stochastic, 
Eq. (24) which governs the transport of the joint scalar PDF is deterministic. 


While Eq. (41) dictates the spatial evolution of the FDF in the physical domain, the com- 
positional makeup of the particle evolves simultaneously due to the actions of mixing and 
reaction. This is accomplished by solution to the differential equation 


dt 


d / d<f> a d(j> a | \ 


(42) 


where the quantity </>+ = <f> Q (Xi(t),t) denotes the scalar value of the particle with the La- 
grangian position vector X t . With the introduction of the LMSE model, Eq. (42) takes on 
the form 

~ <Aa) + (43) 

Effectively, Eqs. (41) and (43) constitute an equivalent system with the same FDF that is 
governed by Eq. (28). It should be clear that while each individual particle evolves in the 
physical domain in a stochastic manner, the statistics of these particles which are used to 
evaluate the FDF evolve in a deterministic manner according to Eq. (28). The essence of 
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the Monte Carlo approach is not to directly simulate fluid particles and the associated FDF, 
but rather to indirectly determine the FDF by tracking notional particles which evolve with 
the identical statistics. 


4 Numerical Approach 

In the numerical implementation of the equivalent system discussed in the previous section, 
the FDF is represented by N p Monte Carlo particles, each with a set of scalars <^v*)(X' n )(f ), t) 
and Lagrangian position vector X^ n E Numerically, a splitting operation is used to treat 
the effects of physical transport, molecular mixing and chemical reaction separately. The 
simplest means of simulating Eq. (41) is via the Euler-Maruyamma approximation (Kloeden 
and Platen, 1995). This preserves the Markovian character of the diffusion process (Gillespie, 
1992; Papoulis, 1965; Billingsly, 1979; Helfand, 1979; Schuss, 1980; Ross, 1996) and facilitates 
affordable computations. Operationally, this implies that if the time interval [to, t > to] is 
divided into increments to, ti, t 2 , ... f,, t,+x,. .. t with to < t\ < t 2 < . . . t; < L+i <...t, the 
“conditional” probability (for an event A) obeys: 

Prob. {M,t n |t 1 ,t 2 ,...t ni } = Prob. {*4, (44) 

Using the Euler-Maruyamma scheme, Eq. (41) takes on the discretized form: 

x?(t k+1 ) = X?(t k ) + DUhW + B n (t k )(Aty' 2 C(h) (45) 

where D™{t k ) = D t (X^(t k )) and B n (t k ) = B(X.^(t k )). Higher order numerical schemes for 
solving Eq. (41) are available (Kloeden and Platen, 1995), but one must be very cautious 
in actual simulations. It must be recognized that since in LES, the diffusion term in Eq. 
(38) strongly depends on the stochastic process X(t), the numerical scheme must preserve 
the Ito-Gikhman nature of the process. The coefficients D t and B require the input of the 
filtered mean velocity and the diffusivity (molecular and subgrid eddy). These quantities are 
determined by conventional LES at standard finite difference grid points. Since the Monte 
Carlo particles are not restricted to the LES grid points, fourth order Lagrange polynomials 
are utilized to interpolate the desired quantities to the particle. 
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At each time step, the compositional values are subject to change due to subgrid molecular 
mixing and chemical reaction. Eq. (43) may be integrated numerically to simulate the 
effects of chemical reaction and dissipation simultaneously. Alternately, Eq. (43) may be 
treated in a time-split manner in which the effect of mixing and chemical reaction are treated 
separately. The advantage to this approach is that the compositional changes due to mixing 
may be expressed in an analytical manner. With the reaction term set to zero, Eq. (43) 
may be integrated analytically over a time step At to determine the change in composition 
due to molecular mixing. Consequentially each particle changes according to 


(«)”" = <K + (« - 


1 


C^mAt 


(46) 


at each time step. After this step, the particles undergo reaction. This is performed readily 
by sweeping over all the particles and determining the fine grain reaction rates and 
modifying the composition of each of the elements accordingly: 


<j>:( t +At) = (<j>:r*+u n a At. 


(47) 


The mixing procedure given by the LMSE model implementation Eq. (46) requires filtered 
value of the compositional vector <f> a as an input. The estimation of subgrid moments at a 
given point is conducted by consideration of particles within some volume centered at the 
point of interest. Effectively, this finite volume constitutes an “ensemble domain” (not to 
be confused with the “filter domain” characterized by the filter length Ag) in which the 
FDF is represented discretely by numerous stochastic particles. This is necessary as with 
probability one no particles will coincide with the point in question. Essentially, each particle 
represents a “fine grain value” and may be considered a single realization of the flow. Strictly 
speaking, it is not correct to literally equate a stochastic particle to a fluid sample, however 
the collective statistics they are used to generate are equivalent. In the present work, a 
simple box average is utilized as the ensemble domain at the finite difference grid nodes, and 
the fourth order Lagrangian interpolation procedure is used to interpolate these means to 
the particle positions. The size of the ensemble domain is an important issue and pertinent 
results are presented in a later section. Other methods used to calculate moments, such as 
those based on cubic splines, are not a subject of the current investigation. The numerical 
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scheme only relies upon the input of the filtered scalar value, and not its derivative, in the 
mixing model implementation. For this purpose the simple volume averaging procedure is 
sufficient. 

The numerical procedure used to calculate the mean hydrodynamic transport variables uti- 
lizes a high order finite difference scheme as discussed by Carpenter (1990). This scheme is 
a variant of the well known MacCormack scheme in which fourth order compact differences 
are used to discretize the spatial derivatives. A second order symmetric predictor-corrector 
sequence is employed to achieve second order accuracy in time. The code solves the Navier 
Stokes equations in fully compressible form, however for the purposes of this paper all sim- 
ulations are conducted at a low Mach number (M — 0.3) so as to minimize the effects of 
compressibility. The nature of the finite difference scheme is independent of the Monte Carlo 
method and alternative methods could be used in its place. One attractive feature of the 
Lagrangian scalar FDF Monte Carlo scheme is that it could be incorporated rather easily in 
the extensive stockpile of existing fluid dynamic codes. 


5 Results 

To demonstrate the effectiveness of the Monte Carlo FDF approach, two flow configurations 
are considered: (1) a temporally developing mixing layer and (2) a planar jet. For both 
configurations, a constant rate, non heat releasing chemical reaction of the type A + B — ► P is 
simulated. In this study, only two dimensional flows are considered. The fundamental theory 
is three dimensional and hence no special treatment is required to extend the simulations to 
3 spatial dimensions. 

The mixing layer configuration consists of two coflowing streams traveling at different ve- 
locities and merging at the trailing edge of a partition plate. Two reactants, A and J9, are 
introduced into the high and the low speed streams, respectively. To expedite the formation 
of large scale vortices, low amplitude perturbations are initially superimposed on the mean 
flow. The flow field which develops in this setting is dominated by large scale coherent struc- 
tures. In the spatial jet, reactant A is injected into the domain in the high speed stream 
while reactant B enters the domain in the lower speed coflow. Disturbances are introduced 
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at the inlet plane so as to facilitate the formation of large scale, energy containing eddies. 

The assessment of the models is facilitated by direct comparison with DNS data. The resolu- 
tion in DNS is dictated by the magnitudes of the physical parameters, with sufficient testing 
on the independency of the results to the number of grid points. The highest resolution 
in temporal shear layer simulations consists of 434 x 577 grid points. With this resolution, 
reliable DNS with a Reynolds number Re — 2800 and a Damkohler number Da = 2 (based 
on the velocity difference and the vorticity thickness at initial time) are possible. The length 
in streamwise direction is chosen to two times the wavelength of the most unstable mode 
given by linear stability theory. The scenario of the simulation showed the rollup of the span- 
wise vorticity, resulting in two spanwise rollers. Subsequently, pairing of these vortices was 
observed, reducing the number of vortices to one in later stages. The DNS of the spatially 
evolving jet was conducted on an evenly spaced grid with resolution of 601 x 301 grid points. 
With this resolution, accurate simulations up to Re = 10,000 and Da = 2 (based upon the 
jet diameter and centerline velocity) are possible. Large scale structures are generated by 
forcing of the cross-stream velocity at the inlet plane at a frequency corresponding to the 
most unstable mode as determined by linear stability theory. The velocity ratio of the high 
to low speed streams is 0.5. The computational domain extends 14 diameters downstream 
and is 7 diameters wide. In both flows, all the species {A,B, P) are assumed thermodynam- 
ically identical and the fluid is assumed to be calorically perfect. The value of the turbulent 
Schmidt number is set equal to 0.7 in all simulations. 

The spatial filter A g is set to twice the grid spacing for all simulations. The effect of the 
size of the “ensemble domain” is considered by comparing the results for FDF simulations 
performed using uniform box domains of dimension Ax x Ay and 2Ax x 2Ay (in the present 
simulations the grid spacing Ax = Ay is constant throughout the computational domain). 

In order to assess the performance of the FDF methodology in predicting the scalar cor- 
relations in the chemical source terms, finite difference LES simulations are conducted in 
which the reaction rate is assumed to be a function of the mean species concentrations, 
thus neglecting any correlation due to the small scales. The convective turbulent transport 
terms are modeled in a “consistent” manner by utilizing the same eddy diffusivity, so any 
discrepancy can be attributed to the neglection of the scalar correlations in the mean reac- 
tion rates. In the degenerate case of an nonreactive scalar where the species conversion rate 
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is zero, the moment equation solved by the finite difference should yield the same results 
as those generated by solution to the FDF transport equation. Any deviation between the 
two must be attributed to differences in the numerical solution (i.e. Monte Carlo vs. finite 
difference). 

5.1 Non-reacting simulations 

In order to assess the validity of the convection, diffusion and dissipation numerical pro- 
cedures, Monte Carlo FDF and “conventional” finite difference LES simulations were per- 
formed for a conserved scalar. As noted above, it is expected that the difference between 
the two approaches should be small, any discrepancy owing only to the different numerical 
methodologies. 

In Fig. 1, results are presented of the LES of the non-reacting temporally developing mixing 
layer. The contour plots of the mean conserved scalar as obtained by (a) the FDF and (b) 
conventional LES approaches illustrate the consistency between the FDF and conventional 
LES methodologies. In fact, the Monte Carlo results are even more attractive due the 
Lagrangian nature of the solution procedure. Note the oscillations in the results provided by 
the conventional finite difference LES approach. While the finite difference method suffers 
from potential over and undershoots (which is controlled somewhat by clipping non physical 
values) in lower resolution LES simulations, there is no numerical mechanism to cause such 
errors in the Monte Carlo scheme (at least in the case of non-reacting scalars). Even with the 
addition of the subgrid diffusivity at this low resolution simulation, the conventional finite 
difference LES scheme is not capable of providing a solution free of numerical noise. 

Figs. 2-3 show the variation of the mean (filtered) value and generalized variance of the 
conserved scalar, <f>, across the shear layer. These statistics are gathered by averaging the 
data over the homogeneous (x) direction. This is represented by [] x . Such statistics in 
the temporal simulations are the counterpart to Reynolds averaged quantities for spatially 
evolving flows. The most significant difference between the results obtained from FDF and 
conventional LES method are at the final time of the simulation ( t « 44) when the merging 
of the two pairs of vortices is completed. 
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Figure 2(a) shows that the values of [<f\ x evaluated from finite difference and Monte Carlo 
methods are indeed very close even when the FDF sample domain size is relatively large 
(2Ax x 2Ay). The results for [<j) 2 ] x (not shown) exhibit a similar behavior. The difference 
between FDF and LES is better observed in Figs. 2(b) and 2(c) where the variations of 
[a 2 ] x = [ <f> 2 — (f> 2 ] x with y for different mixing frequency constant (C^) are shown. As ex- 
pected with increasing C^, the variance is significantly decreased. Nevertheless, the deviation 
between the variances calculated from FDF and LES methods is consistently insignificant. 
The numerical error due to the oscillatory behavior in the finite difference method is one of 
the reasons for the observed deviation between variances. However, the difference is mainly 
due to the finite size of the ensemble domain in the FDF calculation and is decreased as the 
size of ensemble domain for FDF calculation is decreased (Figs. 2(b) and 2(c)). Ideally, the 
best results from the Monte Carlo procedure are obtained when the size of sample domain 
is infinitely small and the number of particles within this domain is infinitely large. Since 
it is necessary to deal with a finite number of particles, a smaller domain size results in a 
fewer particles to construct the FDF with the consequence of higher statistical error. Alter- 
nately, a larger ensemble domain will result in a decrease in statistical error, however the 
consequence is an increase in spatial error which manifests itself in a diffusive type effect. 
The optimum size varies with many factors and can not in general be specified a priori. The 
results in Figs. 2(b) and 2(c) show that the variances evaluated from LES and FDF are 
reasonably well compared when the ensemble domain is Ax x Ay. Therefore in most of 
the simulations this sample size is used. Other results shown below suggest that important 
statistical quantities such as the product thickness are not very sensitive to the sample size 
for calculation of FDF. 

The sensitivity of the results to the number of particles and the time step of calculations are 
shown in Fig. 3. The time step may be especially important in the numerical integration 
of the stochastic Langevin equation and is decreased by decreasing the magnitude of CFL 
number. The value of [cr 2 ] x does not vary significantly with the particle number density as 
long as the particle density is not significantly low. This is observed in Fig. 3(a), where it is 
shown that even at the final time of the simulation [cr 2 ) x values calculated from particles are 
almost invariant with respect to particles number density. It is also shown in Fig. 3(b) that 
by decreasing the time step below that is required by stability considerations [ cr 2 ] x remains 
invariant. The results for [<j>] x (not shown) exhibits much less sensitivity to the initial number 
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of particles per cell and also to the time step. 

FDF and finite difference simulations of the spatially evolving jet for a conserved scalar 
further support the observations noted in the case of the temporal shear layer. Fig. 4. 
are plots of the time averaged conserved scalar as a function of the cross-stream coordinate 
for several downstream positions. It is evident that the mean scalar is predicted in nearly 
identical form for the FDF and finite difference LES approaches. As in the case of the 
temporal shear layer, the results of the generalized variance exhibit a somewhat sensitive 
behavior to the size of the ensemble domain. This is confirmed by Fig. 5 which shows the 
variation of the time averaged generalized subgrid variance with the cross-stream direction 
for 2 downstream positions. The deviation between the FDF and finite difference LES 
approaches diminishes as the ensemble domain decreases. 

5.2 Reactive simulations 

The primary motivation for utilizing the FDF approach is that the scalar correlations appear 
in closed form thus avoiding the closure problem. The results of the previous section illus- 
trate the equivalence of the FDF and moment closure approaches in the non-reactive case. 
This equivalence was exploited to demonstrate the effectiveness of the numerical scheme in 
predicting the combined effects of convection, diffusion and mixing. In this section the effects 
of nonequilibrium chemical reaction are considered. Under these conditions a consistently 
between the FDF and moment closure approaches no longer exists, and the full benefits of 
the FDF approach are be realized. 

To conduct a resolution study, model free simulations at different resolutions are conducted. 
Temporal evolution of the vorticity and product thicknesses corresponding to these cases are 
shown in Fig. 6. In Fig. 6(a) it is observed that the calculated vorticity thickness for the 
simulation with 290 X 385 grid points is identical with that calculated based on 434 x 577 grid 
points. However, by decreasing the resolution to 38 x 49 the vorticity growth is decreased. 
Similar results are shown in Fig. 6(b) for the product thickness. Even at Da=2, the product 
thickness increases only very slightly when the resolution is decreased from 434 x 577 to 
290 x 385. On the other hand, the predicted value of S product at low resolution (38 x 49) 
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is much higher than that obtained from high resolution simulations. Based on the results 
shown in Fig. 6, a resolution of 434 x 577 grid points is sufficient to resolve the hydrodynamic 
and scalar fields. The FDF and finite difference LES simulations are conducted on a 38 x 49 
grid. 

It is evident from the results shown in Fig. 6 that in LES the unresolved subgrid scales have a 
significant effect on the resolved scales. To show this and also to demonstrate the importance 
of tab in LES, in Fig. 7 we consider the variation of the tab (x-averaged) and its subpart, 
Rab = Wb-J'aJ'b along y direction. Following Leonard (1974) and Germano (1986), t A b 
can be decomposed into subparts. Each of the terms in this decomposition is expected to 
constitute a significant part of the SGS unmixedness when a filter allowing overlap between 
resolved and unresolved scales is used, as for instance the Gaussian or the box filter. One 
of these subparts is R A b, the “Reynolds” part. The results in Fig. 7 show that throughout 
the simulation R A b is only a fraction of t A b ■ This suggests that the modeling of t a b in 
LES is more difficult than the modeling of the corresponding term in RANS calculations. 
In the previous LES of reacting flows t A b is either effectively ignored (Schumann, 1989; 
Liou et al., 1994) or is modeled in an ad hoc basis (Fureby and Lofstrom, 1994). 

In Fig. 8, the values of 6 product calculated from FDF method are compared to that obtained 
from filtered DNS data. In this calculation the Smagorinsky model (rather than the MKEV 
model) is used to evaluate the subgrid viscosity. The results in Fig. 8 show that regardless 
of the magnitude of the Smagorinsky constant (C s ), S pToduct is not well predicted by FDF 
method even though the chemistry in this method is closed. The main reason is the inade- 
quacy of the Smagorsinky model to predict the correct value of the subgrid viscosity in this 
transitional flow, which in turn results in a poor prediction of the mixing frequency. While 
the reaction is treated exactly in the FDF approach, the composition at which this reaction is 
evaluated at is effected by the selection of the mixing frequency. This points the importance 
of proper evaluation of the mixing frequency. It is well established that the Smagorinsky 
closure generates excessive damping on the resolved large scales in transitional regions and 
consequently wrong prediction of the growth rate of the shear layer. In the temporal shear 
layer C s should be initially set to zero because the flow is laminar at this time. Gradually 
it should increase with the growth of the layer. As shown in Fig. 8, with choosing the 
Smagorinsky coefficient to be a linear function of time, the transitional nature of the flow is 
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well captured as well as the value of the subgrid viscosity, and in turn the mixing frequency. 
As a result, the values of 6 pro d U ct calculated with C s oc t compare favorably with the DNS 
values. Indeed, it is intriguing how well & pr oduct is predicted with this simple functional form 
for C s . No attempt is made to propose this form of C s for any flow, rather the intention is 
to emphasize the importance of the subgrid viscosity and particularly the associated mixing 
frequency in the FDF method and provides motivation to find a better model to calculate 
the subgrid viscosity. 

Motivated to better predict the eddy viscosity, the MKEV model, Eq. (12), is adopted. The 
improved prediction of the subgrid viscosity has favorable consequences. First, it results 
in better prediction of the hydrodynamics which affects the mixing and reaction processes. 
Second, it yields better predictions of the subgrid diffusivity in the Langevin equation and 
the mixing frequency in the LMSE model. In Fig. 9 it is shown that the temporal evaluation 
of the vorticity thickness calculated from the above MKEV model is closer to that calculated 
from DNS than those calculated from Smagorinsky model. The Smagorinsky model with 
time varying (linear in time) coefficient better compares with the DNS than the Smagorinsky 
model with constant coefficient. Proper prediction of the subgrid viscosity has an even more 
significant effect on the chemistry as a result of the improvement of the mixing frequency. 
This is observed in Fig. 10(a) where S pro d uc t as calculated by the FDF method is shown. 
It is observed that the FDF simulation compares favorably with the DNS over a range of 
Damkohler numbers. In Fig. 10(b) similar results are shown for LES calculations with 
models for the SGS stresses and scalar fluxes and no model for SGS unmixedness. Similar to 
that observed in Fig. 6(b), the LES values of S product are much higher than the corresponding 
DNS values, suggesting that the important effect of SGS unmixedness cannot be overlooked. 

In Fig. 11(a) the variation of the streamwise averaged product across the shear layer eval- 
uated from DNS, LES (with no model for the subgrid unmixedness) and FDF is shown. 
Consistent with the results shown in Fig. 10, it is observed that the FDF method reason- 
ably well predicts the product distribution while the values calculated from LES method 
are significantly higher than the DNS values. The trend is similar at different Damkohler 
numbers. Furthermore it is shown in Fig. 11(b) that if the reaction source term is incorrectly 
calculated to be a function of the mean scalar values (<h a = Dacf> A <j> B ) in the FDF approach, 
the corresponding results for [<j) P ] x are similar to that evaluated from the LES method with 
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no reaction model. This is in accordance with the consistency of the Monte Carlo procedure 
and finite difference method as established above. The slight difference between the product 
profiles obtained from DNS and FDF methods in Fig. 11(a) is mainly due to inadequacy of 
the eddy viscosity type of closures to predict the velocity field. As a result, the transport of 
scalars (including the product) is not well captured by FDF (or LES) method. This issue is 
further discussed below where the possible remedies to resolve it is suggested. 

The reasonably good prediction of the product formation by the FDF approach is due to 
capability of the method to explicitly calculate the nonlinear reaction source term. It is 
important that the reaction rate be correctly predicted by the model throughout the scalar 
evolution. In Fig. 12 the reaction source term calculated from the DNS, LES and FDF 
methods are shown at two different times. Consistent with the results in Fig. 11, the 
reaction rate term obtained from FDF method compares favorably with that of the DNS. 
The predicted reaction rate evaluated from LES method is much higher than the DNS 
values because the SGS unmixedness term is essentially ignored. The corresponding SGS 
unmixedness term is explicitly accounted for in the FDF method. 

The robustness of the FDF method is further established in Fig. 13 where the sensitivity of 
the method to different parameters is tested. The results show that the product evaluated 
from the FDF method is essentially invariant to the number density of the particles provided 
that this density is not very low. Furthermore, it is observed that the chemical reaction is 
insensitive to the size of sample that is used to evaluate FDF. While the mean value of the 
scalar used in the mixing model for a given particle should be the mean value of the scalar 
at the particle location, it was observed that the mean value at the nearest finite difference 
grid point could be substituted. This eliminates the need for interpolating the mean scalar 
field to the particle locations. The result in both cases is essentially the same. 

It was shown in Fig. 2(a) that the mean value of a conserved scalar <j> calculated from FDF 
and LES methods are very close. As established mathematically, this similarity should always 
exist regardless of the values of mixing and flow parameters and also the type of models that 
are used for SGS stresses and scalar fluxes. This is confirmed in Fig. 14(a) where [<f > A ] r for 
Da = 0 is shown. Additionally it is observed in this figure that at initial period of scalar 
evolution the LES and FDF values compare quite well with DNS values. However at latter 
times the difference between DNS and FDF (or LES) results is increased. The main reason 
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for this difference as explained before and as shown in Fig. 14 is the damping of the growth 
of shear layer that is caused by eddy diffusivity closures. Consistently, the reactant scalars 
also show a difference (although not very significant) between FDF and DNS results at later 
times of flow development (Fig. 14(b)). Expectedly, the deviation of the mean reactant 
values calculated from LES method with respect to DNS values is much more than that of 
the FDF. 

The effectiveness of the FDF approach to predict the slightly more complex jet configuration 
is demonstrated in Fig. 15, where contours of the product mass fraction are plotted for the 
FDF, filtered DNS and LES (without reaction model). A resolution 151 x 76 is used in the 
FDF and LES simulations. A total of 20 particles per grid cell (yielding 80 particles per 
ensemble domain of 2Ax X 2A y) were initialized at the start of the calculation. In order 
to reduce the computational expense, stochastic particles are only initialized in the region 
extending from y = —1.75 D to y = 1.75 D where y is the cross-stream coordinate and D 
is the jet diameter. New particles are introduced at the inlet only in this region at a rate 
proportional to the local flow velocity to maintain a constant particle density. It is assumed 
no reaction takes place in cells that do not contain particles. This approach requires that 
particles only be seeded in areas where the flow is undergoing chemical reaction yielding a 
considerable savings. The DNS is better represented in the FDF simulation than in the LES 
simulation, where the product formation in the coherent structures in considerably higher. 
This is further evidenced by Fig. 16 which shows the spatially evolving product thickness 
plotted as a function of the downstream coordinate. Additionally, Fig. 17 exhibits the 
time averaged product profiles as a function of the cross-stream coordinate at a downstream 
position of 2, 7 and 11 jet diameters. A critical assessment of these figures clearly indicates 
that neglect of the SGS unmixedness in the LES procedure results in gross errors in product 
formation while the FDF approach is able to resolve these effects and . better predict the DNS 
data. 
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6 Conclusions 


A filtered density function (FDF) method suitable for chemically reactive flows is developed 
in the context of large eddy simulation. The advantage of the FDF methodology is its in- 
herent ability to resolve SGS scalar correlations that otherwise have to be modeled. Because 
of the lack of robust models to accurately predict these correlations in turbulent reactive 
flows, simulations involving turbulent combustion are often met with a degree of skepticism. 
The FDF methodology avoids the closure problem associated with these terms and treats 
the reaction in an exact manner. The scalar FDF approach is particularly attractive since 
it can be coupled with existing hydrodynamic CFD codes. 

For conserved scalars the consistency between the FDF and moment closure approaches is 
demonstrated. This is observed by integrating the FDF transport equation to derive the 
corresponding moment equations. This consistency no longer exists in the reactive case and 
the benefits of the FDF approach are realized. No model is required in the FDF approach 
and thus the equivalent moment closure does not exist. The principle of equivalent systems 
is utilized in order to develop a tractable solution technique to indirectly calculate the FDF. 
A Lagrangian Monte Carlo method is employed to generate an equivalent system with the 
same FDF as that of true fluid particles. 

Comparison of the FDF solution to that calculated by a conventional finite difference LES 
simulation utilizing a consistent moment closure yields favorable results for a conserved 
scalar. Such comparison is necessary to build confidence that the numerical particle method 
is capable to accurately capture the effects of convection, diffusion and mixing. Simulations 
of chemically reactive flows demonstrate the effectiveness of the FDF method in capturing 
the mean reaction rate. Comparison with a finite difference solution which does not attempt 
to model the SGS scalar covariance indicates the neglect of this term leads to unphysically 
high product conversion rates. Comparison of results between the FDF approach and those 
generated by DNS of two dimensional reacting shear layers and jets demonstrate that the 
FDF approach is capable of accurately capturing the primary features of unsteady combus- 
tion. 

Although the present methodology has been developed for isothermal reactions, the extension 
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to exothermic variable density flows imposes no serious obstacles and is currently under 
investigation by the authors. With the generalization to reactive flows with heat release, 
arbitrarily complicated reaction schemes may be implemented with little difficulty since no 
special treatment is necessary for the source terms in the FDF approach. The scalar domain 
must be modified to include the temperature (or energy) so that any correlations involving 
temperature may be evaluated. Utilizing such a scalar FDF approach, it is conceivable that 
LES of complex reactive flows with realistic chemical kinetics may be routinely simulated 
for engineering applications in the near future. The fact that the scalar FDF approach 
avoids the closure problem for the source terms is significant as no model has yet been 
proposed for use in LES which has universal applicability. This is the primary stumbling 
block for implementation of LES to evaluate unsteady combustion processes. Even for the 
simple chemistry scheme considered in this study there are many difficulties in evaluating 
the nonlinear reaction source term by conventional LES methods. This term is closed in the 
FDF approach and imposes no difficulty even if the reaction takes on a more complex form. 

One drawback of the scalar FDF approach is that any correlations involving the velocity 
field (such as SGS stresses and SGS mass and heat fluxes) must be modeled using existing 
gradient diffusion models. This drawback may be readily overcome by considering the joint 
velocity-scalar FDF in which the velocity correlations appear in closed form eliminating the 
dependence on such models (Pope, 1985). This area has seen considerable development in 
the area of RANS (Pope, 1994). Improvement in the prediction of the subgrid viscosity is 
also necessary (and perhaps more significant) in the estimation of the turbulent frequency in 
the molecular mixing model. The extension of FDF method to include the velocity domain 
is also the subject of current investigation. 
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Fig. IS: Product mass fraction contours (a) DNS 

(b) FDF and (c) LES with no reaction model . 

















